#III级杂音特征图
from python_speech_features import fbank
from python_speech_features import mfcc
from python_speech_features import logfbank
import os
import shutil
import re
import wave
import librosa
import matplotlib.pyplot as plt
import librosa.display
import soundfile
import numpy as np
def FrameTimeC(frameNum, frameLen, inc, fs):
    ll = np.array([i for i in range(frameNum)])
    return ((ll - 1) * inc + frameLen / 2) / fs
def FrequencyScale(nfilt,fs):
    high_freq_mel = (2595 * np.log10(1 + (fs / 2) / 700))  # 求最高hz频率对应的mel频率
# 我们要做40个滤波器组，为此需要42个点，这意味着在们需要low_freq_mel和high_freq_mel之间线性间隔40个点
    mel_points = np.linspace(0, high_freq_mel, nfilt + 2)  # 在mel频率上均分成42个点
    hz_points = (700 * (10 ** (mel_points / 2595) - 1))  # 将mel频率再转到hz频率
    return hz_points
wav_file='E:/HZH/heart_data/2530_AV.wav'
filename="clinic_data/wav_cut_data/0081_PV_Soft_1.wav"
filename1="clinic_data/wav_cut_data/0134_Erb_Loud_1.wav"
filename2="clinic_data/wav_cut_data/0092_PV_Absent_2.wav"
filename3="clinic_data/wav_cut_data/x0017_PV_Absent_1.wav"
x, fs = librosa.load(filename3,sr=4000)
# x = librosa.resample(x, orig_sr=fs, target_sr=4000)fs=4000
fbank_feat=logfbank(x,fs,winlen=0.025,winstep=0.0125,nfilt=32,nfft=512,lowfreq=0,highfreq=800)#计算对数fbank特征（）
# mfcc_feat=mfcc(x,fs,winlen=0.025,winstep=0.0125,numcep=13,nfilt=32,nfft=512)#计算mfcc特征，（倒谱数，帧数）
print(fbank_feat.shape)
# print(mfcc_feat.shape)
wlen=fs*0.025
feat=fbank_feat.T
print(feat.shape)
freq=FrequencyScale(30,1600)
print("len(freq)",len(freq))
# print("freq:",freq)
frameTime = FrameTimeC(feat.shape[1], wlen, wlen*0.5, fs)
print(len(frameTime))
plt.pcolormesh(frameTime, freq,feat)
plt.xlabel('time(s)')
plt.ylabel('frequency(Hz)')
plt.colorbar()
#plt.savefig('语谱图22.png')
plt.show()